1 Scoping Review

This document visualizes the data for the scoping review.

1.1 Study Designs Across Domains

This plot shows the number of studies using a longitudinal or serial cross-sectional design across several domains.s

### plot domains with study design
### Function to create a single horizontal bar plot for counts of 1
create_combined_horizontal_bar <- function(df) {
  columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
  
  # Create a long-format data frame for ggplot
 df_long <- df %>%
  tidyr::gather(key = "variable", value = "value", columns_of_interest) %>%
  filter(value == 1)
  
  # Arrange levels of 'variable' by count in ascending order
  df_long$variable <- factor(df_long$variable, levels = names(sort(table(df_long$variable), decreasing = FALSE)))
  
  # Define custom colors using hex color codes
  custom_colors <- c("#440154", "#22A884")  # Hex color codes without alpha
  
  # Plotting the horizontal bar plot using ggplot2
  ggplot(df_long, aes(x = value, y = variable, fill = study_design)) +
    geom_bar(stat = "identity", position = position_stack(reverse = TRUE), color = NA) +  # Reverse the stacking order
    scale_fill_manual(values = setNames(custom_colors, unique(df_long$study_design))) +  # Apply custom colors to study_design levels
    theme_minimal() +
    theme(panel.grid = element_blank(), panel.border = element_blank()) +  # Remove grid lines and borders
    labs(x = "Number of Studies", y = NULL, title = NULL) +  # Remove y-axis label and title
    xlab("Number of Studies")  # Set x-axis label
}
# Apply the function to df_final
create_combined_horizontal_bar(df_final)

1.2 Measuring Risk Perception Across Domains

This plot shows the number of studies and how many times they measured risk perception across several domains.

###plot domain with measurement points-----
#create measurement category with 2, 3, and 3+
df_final <- df_final %>%
  mutate(measure = ifelse(times_measured_1 == 2, "2",
                          ifelse(times_measured_1 == 3, "3", "3+")))

# Function to create a single horizontal bar plot for counts of 1
create_combined_horizontal_bar <- function(df) {
  columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
  
  # Create a long-format data frame for ggplot
  df_long <- df %>%
    gather(key = "variable", value = "value", columns_of_interest) %>%
    filter(value == 1)  # Filter only rows where the value is 1
  
  # Arrange levels of 'variable' by count in ascending order
  df_long$variable <- factor(df_long$variable, levels = names(sort(table(df_long$variable), decreasing = FALSE)))
  
  # Plotting the horizontal bar plot using ggplot2 with Viridis color scale
  ggplot(df_long, aes(x = value, y = variable, fill = measure)) +
    geom_bar(stat = "identity", position = position_stack(reverse = TRUE), color = NA) +  # Reverse the stacking order
    scale_fill_viridis_d() +  # Use Viridis color palette for measure variable
    theme_minimal() +
    theme(panel.grid = element_blank(), panel.border = element_blank()) +  # Remove grid lines and borders
    labs(x = "Number of Studies", y = NULL, title = NULL) +  # Remove y-axis label and title
    xlab("Number of Studies")  # Set x-axis label
}

# Apply the function to df_final
create_combined_horizontal_bar(df_final)

1.3 Risk Formulation Across Domains

This plot illustrates the different risk formulations used to measure risk perception across various domains.

#kind of risk (risk, concern, worry)
create_combined_horizontal_bar <- function(df) {
  columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
  
  # Create a long-format data frame for ggplot
  df_long <- df %>%
    gather(key = "variable", value = "value", columns_of_interest) %>%
    filter(value == 1)  # Filter only rows where the value is 1
  
  # Plotting the horizontal bar plot using ggplot2 with Viridis color scale
  ggplot(df_long, aes(x = variable, fill = as.factor(risk_1))) +
    geom_bar(position = "fill", stat = "count", color = "white") +
    coord_flip() +  # Flip the coordinates for a horizontal bar plot
    theme_minimal() +
    theme(panel.grid = element_blank(), panel.border = element_blank(),
          axis.text.x = element_blank(), axis.title.x = element_blank()) +  # Remove x-axis labels and title
    labs(x = NULL, y = NULL, title = NULL) +  # Remove axis labels and change title
    scale_fill_viridis_d(name = "Risk Formulation")  # Use Viridis color palette and customize legend title
}

# Apply the function to df_final
create_combined_horizontal_bar(df_final)

1.4 Event Occurence Across Domains

This plot shows the occurrence of events across various domains.

# Define the function
create_combined_horizontal_bar <- function(df) {
  columns_of_interest <- c('health', 'finance', 'political', 'crime', 'nature', 'nuclear', 'social')
  
  # Create a long-format data frame for ggplot
  df_long <- df %>%
    gather(key = "variable", value = "value", columns_of_interest) %>%
    filter(value == 1)  # Filter only rows where the value is 1
  
  # Create a new factor combining intervention and exposure columns into "event" and "no event"
  df_long$event_category <- ifelse(
    df_long$intervention_yesno_1 == 0 & df_long$exposure_yesno_1 == 0, 
    "No Event", 
    "Event"
  )
  
  # Plotting the horizontal bar plot using ggplot2 with Viridis color scale
  ggplot(df_long, aes(x = variable, fill = event_category)) +
    geom_bar(position = "fill", stat = "count", color = "white") +
    coord_flip() +  # Flip the coordinates for a horizontal bar plot
    theme_minimal() +
    theme(panel.grid = element_blank(), panel.border = element_blank(),
          axis.text.x = element_blank(), axis.title.x = element_blank()) +  # Remove x-axis labels and title
    labs(x = NULL, y = NULL, title = NULL) +  # Remove axis labels and change title
    scale_fill_viridis_d(name = "Event Occurence", labels = c("No Event", "Event"))  # Use Viridis color palette and customize legend title
}

# Apply the function to df_final
create_combined_horizontal_bar(df_final)

1.5 Demographics Table Placeholder

Participant characteristics and demographics.

#Table

2 Temporal Stability

2.1 Model Fitting and Evaluation - Meta-Analytical Stability and Change Model (MASC)

2.1.1 Model 1: Fitting with Correlation, Interval, and Size

# Define family
family <- brmsfamily(
  family = "student",
  link = "identity"
)

# Define the formula
formula_m1 <- bf(
  cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
  nlf(rel ~ inv_logit(logitrel)),
  nlf(change ~ inv_logit(logitchange)),
  nlf(stabch ~ inv_logit(logitstabch)),
  logitrel ~ 1,
  logitchange ~ 1,
  logitstabch ~ 1,
  nl = TRUE
)

# Define the weakly informative priors
priors <- 
  prior(normal(0,1), nlpar="logitrel", class = "b") +
  prior(normal(0,1), nlpar="logitchange", class = "b") +
  prior(normal(0,1), nlpar="logitstabch", class = "b") +
  prior(cauchy(0,1), class = "sigma")

# Fit the model
fit_masc_m1 <- brm(
  formula = formula_m1,
  prior = priors,
  family = family,
  data = df,
  cores = 2,
  chains = 2,
  iter = 6000,
  warmup = 2000,
  control = list(max_treedepth = 10, adapt_delta = 0.95),
  seed = 1299
)
## Compiling Stan program...
## Start sampling
# Model summary
summary(fit_masc_m1)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1) 
##          rel ~ inv_logit(logitrel)
##          change ~ inv_logit(logitchange)
##          stabch ~ inv_logit(logitstabch)
##          logitrel ~ 1
##          logitchange ~ 1
##          logitstabch ~ 1
##    Data: df (Number of observations: 101) 
##   Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Regression Coefficients:
##                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept        0.67      0.17     0.38     1.04 1.00     3867
## logitchange_Intercept    -0.76      0.58    -1.84     0.58 1.00     3667
## logitstabch_Intercept    -0.36      1.09    -2.44     1.75 1.00     3479
##                       Tail_ESS
## logitrel_Intercept        4310
## logitchange_Intercept     2581
## logitstabch_Intercept     3832
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.19      0.02     0.15     0.23 1.00     4261     4046
## nu       17.66     12.01     4.20    50.21 1.00     4216     4488
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Plot conditional effects
plot(conditional_effects(fit_masc_m1), points = TRUE)

# Trace plots and parameter estimates
plot(fit_masc_m1, N = 5, ask = TRUE)

# Posterior predictive checks
pp_check(fit_masc_m1, type = "dens_overlay", ndraws = 100)

pp_check(fit_masc_m1, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)

# LOO and diagnostics
model_loo_m1 <- loo(fit_masc_m1, save_psis = TRUE, cores = 2)
plot(model_loo_m1, diagnostic = "k")

plot(model_loo_m1, diagnostic = "n_eff")

2.1.2 Model 2a: Fitting with Demographics - Age, Age^2, Gender

# Define family
family <- brmsfamily(
  family = "student",
  link = "identity"
)

# Define the formula
formula_m2a <- bf(
  cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
  nlf(rel ~ inv_logit(logitrel)),
  nlf(change ~ inv_logit(logitchange)),
  nlf(stabch ~ inv_logit(logitstabch)),
  logitrel ~ 1 + age_dec_c + age_dec_c2  + female_c,
  logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c,
  logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c,
  nl = TRUE
)

# Define the weakly informative priors
priors <- 
  prior(normal(0,1), nlpar="logitrel", class = "b") +
  prior(normal(0,1), nlpar="logitchange", class = "b") +
  prior(normal(0,1), nlpar="logitstabch", class = "b") +
  prior(cauchy(0,1), class = "sigma")

# Fit the model
fit_masc_m2a <- brm(
  formula = formula_m2a,
  prior = priors,
  family = family,
  data = df,
  cores = 2,
  chains = 2,
  iter = 6000,
  warmup = 2000,
  control = list(max_treedepth = 10, adapt_delta = 0.95),
  seed = 1299
)
## Compiling Stan program...
## Start sampling
# Model summary
summary(fit_masc_m2a)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1) 
##          rel ~ inv_logit(logitrel)
##          change ~ inv_logit(logitchange)
##          stabch ~ inv_logit(logitstabch)
##          logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c
##          logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c
##          logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c
##    Data: df (Number of observations: 101) 
##   Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Regression Coefficients:
##                        Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept         0.49      0.17     0.16     0.85 1.00     5441
## logitrel_age_dec_c        -0.00      0.26    -0.47     0.60 1.00     2283
## logitrel_age_dec_c2        0.22      0.24    -0.00     0.90 1.00     1632
## logitrel_female_c          0.49      0.44    -0.36     1.35 1.00     7709
## logitchange_Intercept     -0.74      0.80    -2.29     0.86 1.00     3341
## logitchange_age_dec_c     -0.33      0.59    -1.65     0.92 1.00     2786
## logitchange_age_dec_c2    -0.03      0.33    -0.61     0.71 1.00     2401
## logitchange_female_c       0.16      0.80    -1.43     1.70 1.00     6218
## logitstabch_Intercept      0.60      1.18    -1.94     2.62 1.00     4362
## logitstabch_age_dec_c      0.16      0.86    -1.57     1.82 1.00     5099
## logitstabch_age_dec_c2    -0.77      0.65    -2.00     0.59 1.00     2338
## logitstabch_female_c       0.02      0.99    -1.94     1.91 1.00     8253
##                        Tail_ESS
## logitrel_Intercept         4035
## logitrel_age_dec_c         2287
## logitrel_age_dec_c2        1000
## logitrel_female_c          5304
## logitchange_Intercept      4888
## logitchange_age_dec_c      2657
## logitchange_age_dec_c2     2167
## logitchange_female_c       4874
## logitstabch_Intercept      5590
## logitstabch_age_dec_c      5164
## logitstabch_age_dec_c2     2477
## logitstabch_female_c       5755
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.18      0.02     0.13     0.23 1.00     3366     3896
## nu       13.96     11.18     2.85    44.29 1.00     3860     4131
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Plot conditional effects
plot(conditional_effects(fit_masc_m2a), points = TRUE)

# Trace plots and parameter estimates
plot(fit_masc_m2a, N = 5, ask = TRUE)

# Posterior predictive checks
pp_check(fit_masc_m2a, type = "dens_overlay", ndraws = 100)

pp_check(fit_masc_m2a, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)

# Leave-One-Out (LOO) and diagnostics
model_loo_m2a <- loo(fit_masc_m2a, save_psis = TRUE, cores = 2)
plot(model_loo_m2a, diagnostic = "k")

plot(model_loo_m2a, diagnostic = "n_eff")

2.1.3 Model 2b: Fitting with Measurement Characteristics - Domain, Item, Event

# Define family
family <- brmsfamily(
  family = "student",
  link = "identity"
)

# Define the formula
formula_m2b <- bf(
  cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
  nlf(rel ~ inv_logit(logitrel)),
  nlf(change ~ inv_logit(logitchange)),
  nlf(stabch ~ inv_logit(logitstabch)),
  logitrel ~ 1 + domain + item + event,
  logitchange ~ 1 + domain + event,
  logitstabch ~ 1 + domain + event,
  nl = TRUE
)

# Define the weakly informative priors
priors <- 
  prior(normal(0,1), nlpar="logitrel", class = "b") +
  prior(normal(0,1), nlpar="logitchange", class = "b") +
  prior(normal(0,1), nlpar="logitstabch", class = "b") +
  prior(cauchy(0,1), class = "sigma")

# Fit the model
fit_masc_m2b <- brm(
  formula = formula_m2b,
  prior = priors,
  family = family,
  data = df,
  cores = 2,
  chains = 2,
  iter = 6000,
  warmup = 2000,
  control = list(max_treedepth = 10, adapt_delta = 0.95),
  seed = 1299
)
## Compiling Stan program...
## Start sampling
# Model summary
summary(fit_masc_m2b)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1) 
##          rel ~ inv_logit(logitrel)
##          change ~ inv_logit(logitchange)
##          stabch ~ inv_logit(logitstabch)
##          logitrel ~ 1 + domain + item + event
##          logitchange ~ 1 + domain + event
##          logitstabch ~ 1 + domain + event
##    Data: df (Number of observations: 101) 
##   Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Regression Coefficients:
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## logitrel_Intercept          1.22      0.37     0.61     2.06 1.00     3725
## logitrel_domainother       -0.14      0.14    -0.40     0.13 1.00     5139
## logitrel_itemsingle         0.43      0.32    -0.07     1.19 1.00     4028
## logitrel_eventevent        -0.68      0.16    -1.01    -0.40 1.00     5017
## logitchange_Intercept      -0.73      0.76    -2.12     0.83 1.00     3907
## logitchange_domainother    -0.43      0.76    -1.83     1.22 1.00     3093
## logitchange_eventevent     -0.31      0.65    -1.58     1.26 1.00     2943
## logitstabch_Intercept       0.04      0.93    -1.84     1.81 1.00     5493
## logitstabch_domainother     0.63      1.00    -1.39     2.39 1.00     3501
## logitstabch_eventevent      0.26      0.95    -1.61     2.10 1.00     3289
##                         Tail_ESS
## logitrel_Intercept          2885
## logitrel_domainother        4837
## logitrel_itemsingle         3159
## logitrel_eventevent         4704
## logitchange_Intercept       5090
## logitchange_domainother     3986
## logitchange_eventevent      2863
## logitstabch_Intercept       5217
## logitstabch_domainother     4837
## logitstabch_eventevent      5142
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.14      0.02     0.09     0.19 1.00     3524     4407
## nu       10.25      9.67     2.26    36.80 1.00     3587     4955
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Plot conditional effects
plot(conditional_effects(fit_masc_m2b), points = TRUE)

# Trace plots and parameter estimates
plot(fit_masc_m2b, N = 5, ask = TRUE)

# Posterior predictive checks
pp_check(fit_masc_m2b, type = "dens_overlay", ndraws = 100)

pp_check(fit_masc_m2b, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)

# Leave-One-Out (LOO) and diagnostics
model_loo_m2b <- loo(fit_masc_m2b, save_psis = TRUE, cores = 2)
plot(model_loo_m2b, diagnostic = "k")

plot(model_loo_m2b, diagnostic = "n_eff")

2.1.4 Model 3: Fitting with Demographics and Measurement Characteristics - Domain, Item, Event, Formulation

# Define the family
family <- brmsfamily(
  family = "student",
  link = "identity"
)

# Define the formula
formula_m3 <- bf(
  cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1),
  nlf(rel ~ inv_logit(logitrel)),
  nlf(change ~ inv_logit(logitchange)),
  nlf(stabch ~ inv_logit(logitstabch)),
  logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
  logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
  logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation,
  nl = TRUE
)

# Define the weakly informative priors
priors <- 
  prior(normal(0, 1), nlpar = "logitrel", class = "b") +
  prior(normal(0, 1), nlpar = "logitchange", class = "b") +
  prior(normal(0, 1), nlpar = "logitstabch", class = "b") +
  prior(cauchy(0, 1), class = "sigma")

# Fit the model
fit_masc_m3 <- brm(
  formula = formula_m3,
  prior = priors,
  family = family,
  data = df,
  cores = 2,
  chains = 2,
  iter = 6000,
  warmup = 2000,
  control = list(max_treedepth = 10, adapt_delta = 0.95),
  seed = 1299
)
## Compiling Stan program...
## Start sampling
# MODEL 3 EVAL: MCMC DIAGNOSTICS --------------------------------------------------------

# Model summary
summary(fit_masc_m3)
##  Family: student 
##   Links: mu = identity; sigma = identity; nu = identity 
## Formula: cor_val | resp_se(sei, sigma = TRUE) ~ rel * (change * ((stabch^interval_val) - 1) + 1) 
##          rel ~ inv_logit(logitrel)
##          change ~ inv_logit(logitchange)
##          stabch ~ inv_logit(logitstabch)
##          logitrel ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
##          logitchange ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
##          logitstabch ~ 1 + age_dec_c + age_dec_c2 + female_c + domain + item + event + formulation
##    Data: df (Number of observations: 101) 
##   Draws: 2 chains, each with iter = 6000; warmup = 2000; thin = 1;
##          total post-warmup draws = 8000
## 
## Regression Coefficients:
##                                Estimate Est.Error l-95% CI u-95% CI Rhat
## logitrel_Intercept                 1.27      0.36     0.67     2.09 1.00
## logitrel_age_dec_c                 0.06      0.15    -0.18     0.41 1.00
## logitrel_age_dec_c2                0.01      0.08    -0.11     0.19 1.00
## logitrel_female_c                  0.24      0.42    -0.57     1.06 1.00
## logitrel_domainother              -0.01      0.15    -0.27     0.30 1.00
## logitrel_itemsingle                0.17      0.29    -0.33     0.81 1.00
## logitrel_eventevent               -0.42      0.15    -0.74    -0.14 1.00
## logitrel_formulationconcern        0.13      0.63    -1.01     1.50 1.00
## logitrel_formulationthreat        -0.32      0.35    -0.99     0.40 1.00
## logitrel_formulationworry          0.80      0.37     0.11     1.59 1.00
## logitchange_Intercept             -0.28      0.89    -1.96     1.50 1.00
## logitchange_age_dec_c              0.38      0.68    -0.92     1.80 1.00
## logitchange_age_dec_c2             0.49      0.61    -0.68     1.79 1.00
## logitchange_female_c               0.06      0.94    -1.85     1.90 1.00
## logitchange_domainother           -0.08      0.78    -1.66     1.47 1.00
## logitchange_itemsingle             0.06      0.95    -1.84     1.85 1.00
## logitchange_eventevent             0.24      0.76    -1.35     1.70 1.00
## logitchange_formulationconcern     0.25      0.94    -1.58     2.13 1.00
## logitchange_formulationthreat      0.57      0.84    -1.16     2.20 1.00
## logitchange_formulationworry       0.30      0.79    -1.24     1.90 1.00
## logitstabch_Intercept              0.04      0.89    -1.73     1.77 1.00
## logitstabch_age_dec_c              0.16      0.61    -1.31     1.19 1.00
## logitstabch_age_dec_c2             0.02      0.47    -1.14     0.73 1.00
## logitstabch_female_c              -0.13      0.91    -1.92     1.73 1.00
## logitstabch_domainother            0.14      0.79    -1.42     1.68 1.00
## logitstabch_itemsingle             0.53      0.91    -1.41     2.18 1.00
## logitstabch_eventevent             0.25      0.78    -1.31     1.74 1.00
## logitstabch_formulationconcern    -0.22      0.90    -1.98     1.53 1.00
## logitstabch_formulationthreat     -0.59      0.89    -2.31     1.22 1.00
## logitstabch_formulationworry      -0.52      0.77    -2.05     1.03 1.00
##                                Bulk_ESS Tail_ESS
## logitrel_Intercept                 3095     3052
## logitrel_age_dec_c                 2569     1699
## logitrel_age_dec_c2                1977     1351
## logitrel_female_c                  6872     5173
## logitrel_domainother               3286     4718
## logitrel_itemsingle                2673     2344
## logitrel_eventevent                3268     4174
## logitrel_formulationconcern        3785     4171
## logitrel_formulationthreat         4152     4037
## logitrel_formulationworry          2849     4353
## logitchange_Intercept              4108     5699
## logitchange_age_dec_c              2066     4377
## logitchange_age_dec_c2             1266     2138
## logitchange_female_c               8640     5622
## logitchange_domainother            3816     4717
## logitchange_itemsingle             2177     4058
## logitchange_eventevent             2106     3880
## logitchange_formulationconcern     7220     5937
## logitchange_formulationthreat      6172     5344
## logitchange_formulationworry       4792     4735
## logitstabch_Intercept              3434     4990
## logitstabch_age_dec_c              1762     3091
## logitstabch_age_dec_c2             1165     1901
## logitstabch_female_c               8073     6119
## logitstabch_domainother            4316     5266
## logitstabch_itemsingle             1904     3466
## logitstabch_eventevent             2048     3512
## logitstabch_formulationconcern     5528     5990
## logitstabch_formulationthreat      5513     5518
## logitstabch_formulationworry       3136     4218
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     0.10      0.03     0.05     0.16 1.00     1675     3489
## nu        4.76      5.30     1.40    19.48 1.00     2101     3951
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# Plot conditional effects
plot(conditional_effects(fit_masc_m3), points = T)

# Trace plots & parameter estimates
plot(fit_masc_m3, N = 5, ask = TRUE)

# MODEL 3 EVAL: PP CHECKS --------------------------------------------------------
pp_check(fit_masc_m3, type = "dens_overlay", ndraws = 100)

pp_check(fit_masc_m3, type = "stat", stat = "mean", ndraws = 1000, binwidth = .001)

pp_check(fit_masc_m3, type = "stat", stat = "sd", ndraws = 500, binwidth = .001)

pp_check(fit_masc_m3, type = "stat", stat = "median", ndraws = 500, binwidth = .001)

pp_check(fit_masc_m3, type = "stat", stat = "mad", ndraws = 500, binwidth = .001)

pp_check(fit_masc_m3, type = "stat_2d")
## Using all posterior draws for ppc type 'stat_2d' by default.

pp_check(fit_masc_m3, type = "scatter_avg")
## Using all posterior draws for ppc type 'scatter_avg' by default.

# MODEL 3 EVAL: LOO --------------------------------------------------------

# LOO & Pareto K
model_loo_m3 <- loo(fit_masc_m3, save_psis = TRUE, cores = 2)
plot(model_loo_m3, diagnostic = "k")

plot(model_loo_m3, diagnostic = "n_eff")

# LOO PIT
w <- weights(model_loo_m3$psis_object)
ppc_loo_pit_overlay(
  y = fit_masc_m3$data$cor_val,
  yrep = posterior_predict(fit_masc_m3),
  lw = w
)
## NOTE: The kernel density estimate assumes continuous observations and is not optimal for discrete observations.

ppc_loo_pit_qq(
  y = fit_masc_m3$data$cor_val,
  yrep = posterior_predict(fit_masc_m3),
  lw = w
)

# MODEL EVAL: LOO - COMPARISON OF ALL MODELS ----------------------------------------------------
loo1 <- loo(fit_masc_m1)
loo2a <- loo(fit_masc_m2a)
loo2b <- loo(fit_masc_m2b)
loo3 <- loo(fit_masc_m3)

loo_compare(loo1, loo2a, loo2b, loo3)
##              elpd_diff se_diff
## fit_masc_m3    0.0       0.0  
## fit_masc_m2b  -7.1       5.3  
## fit_masc_m1  -21.3       6.1  
## fit_masc_m2a -21.7       5.9

2.2 Predictions Time

2.2.1 Plot correlation and interval

nd <- crossing(domain= NA,  
               sei = 0.1,
               item= NA, 
               event=NA,
               formulation=NA,
               age_dec_c= 0,
               age_dec_c2= 0,
               female_c= 0, 
               interval_val=seq(0,5, by = .1))

epred_draws_df <- nd %>% 
  add_epred_draws(fit_masc_m3, re_formula = NA)

ggplot(epred_draws_df) +
  stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) + 
  geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
  theme_minimal() +
  labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
       title = "Behaviour") +
  theme(strip.placement = "outside",
        legend.position = "none",
        legend.margin = margin(-.5,0,0,0, unit="cm"),
        legend.spacing.y = unit(0.15, 'cm'),
        legend.key.width = unit(1, "cm"),
        legend.key.size = unit(.3, "cm"),
        legend.text = element_text(size = 9, color = "grey20"),
        text = element_text(size = 9, color = "grey40"),
        axis.text.y = element_text(vjust=seq(0,1, length.out = 5)),
        axis.text.x = element_text(hjust=c(0,1)),
        title = element_text(size = 9, color = "grey20"),
        plot.tag  = element_text(size = 11, face = "bold", color = "grey20"),
        panel.spacing = unit(.5, "lines"),
        plot.title = element_blank(),
        panel.grid = element_blank(),
        plot.title.position = "plot",
        plot.margin = margin(b = 5, r = 5, l = 5),
        panel.background = element_rect(color = "grey75", fill = NA, size = .4)) +
  scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a")) +
  guides(color = guide_legend(override.aes = list(size = .75)),
         fill = guide_legend(override.aes = list(size = .75)),
         size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0,5)) +
  scale_y_continuous(breaks = seq(0,1,0.25)) +
  scale_x_continuous(breaks = c(0,1,2,3,4,5))

2.2.2 Plot correlation and interval by event

nd <- crossing(domain= NA,  
               sei = 0.1,
               item= NA, 
               event= unique(df$event),
               formulation= NA,
               age_dec_c= 0,
               age_dec_c2= 0,
               female_c= 0, 
               interval_val=seq(0,5, by = .1))


epred_draws_df <- nd %>% 
  add_epred_draws(fit_masc_m3, re_formula = NA)

# epred_draws_dom <- epred_draws_df %>%
#   group_by(interval_val, health_subdomain) %>%
#   mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
#   pivot_wider(names_from = .width, values_from = c(.lower,.upper))


# Set order and capitalise first letter
epred_draws_df$event <- str_to_title(epred_draws_df$event)

epred_draws_df$event<- factor(epred_draws_df$event, levels = c("None","Event"))

df_plot <- df %>% 
  mutate(event = str_to_title(event)) 

df_plot$event<- factor(df_plot$event, levels = c("None","Event"))

# Plot
p_event <- ggplot(epred_draws_df) +
  stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) + 
  geom_point(data= df_plot, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
  # geom_line(data = epred_draws_agg,
  #           aes(x = time_diff_dec*10, y = .epred),
  #           color = "grey95",
  #           size = .5) +
  # geom_line(data = epred_draws_dom, 
  #           aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
  #           color = "#e07f00",
  #           linewidth = .25) +
  # geom_text_repel(data = lbl_dot_df, 
  #                 aes(x = time_diff_dec*10, y = .epred, label = label),
  #                 family = "Source Sans 3", size = 2.5,
  #                 min.segment.length = 0,
  #                 segment.color = "grey50",
  #                 segment.size = .25,
  #                 box.padding = 0.5,
  #                 nudge_x = .5,
  #                 nudge_y = c(.1, -.05)
  # ) +
  facet_wrap(.~event) +
  theme_minimal() +
  labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
       title = "Behaviour") +
  theme(strip.placement = "outside",
        legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
        legend.margin = margin(-.5,0,0,0, unit="cm"),
        legend.spacing.y = unit(0.15, 'cm'),
        legend.key.width = unit(1, "cm"),
        legend.key.size = unit(.3, "cm"),
        # plot.tag.position = c(0,.8),
        legend.text = element_text(size = 10, color = "grey20"),
        text = element_text(size = 12, color = "grey40"),
        axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
        axis.text.x = element_text( hjust=c(0,1)),
        title = element_text(size = 12, color = "grey20"),
        plot.tag  = element_text( size = 12, face = "bold", color = "grey20"),
        panel.spacing = unit(.5, "lines"),
        plot.title = element_blank(),
        panel.grid = element_blank(),
        plot.title.position = "plot",
        plot.margin = margin(b = 5, r = 5, l = 5),
        panel.background = element_rect(color = "grey75", fill = NA, size = .4),
        strip.text = element_text(face = "bold"),
        strip.text.y = element_text(size = 12),
        strip.text.x = element_text(size = 12)) +
  scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a" ))+
  guides(color = guide_legend(override.aes = list(size = .75)),
         fill =  guide_legend(override.aes = list(size = .75)),
         size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
  scale_y_continuous(breaks = seq(0,1,0.25)) +
  scale_x_continuous(breaks = c(0,1,2,3,4,5))


p_event

ggsave(plot = p_event,
       filename = here::here("analysis", "masc_pred_event.jpeg"), 
       dpi = 300, width = 28, height = 30, units = "cm")

2.2.3 Plot correlation and interval by domain

nd <- crossing(domain= unique(df$domain),  
               sei = 0.1,
               item= NA, 
               event= NA,
               formulation= NA,
               age_dec_c= 0,
               age_dec_c2= 0,
               female_c= 0, 
               interval_val=seq(0,5, by = .1))


epred_draws_df <- nd %>% 
  add_epred_draws(fit_masc_m3, re_formula = NA)

# epred_draws_dom <- epred_draws_df %>%
#   group_by(interval_val, health_subdomain) %>%
#   mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
#   pivot_wider(names_from = .width, values_from = c(.lower,.upper))


ggplot(epred_draws_df) +
  stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) + 
  geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
  # geom_line(data = epred_draws_agg,
  #           aes(x = time_diff_dec*10, y = .epred),
  #           color = "grey95",
  #           size = .5) +
  # geom_line(data = epred_draws_dom, 
  #           aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
  #           color = "#e07f00",
  #           linewidth = .25) +
  # geom_text_repel(data = lbl_dot_df, 
  #                 aes(x = time_diff_dec*10, y = .epred, label = label),
  #                 family = "Source Sans 3", size = 2.5,
  #                 min.segment.length = 0,
  #                 segment.color = "grey50",
  #                 segment.size = .25,
  #                 box.padding = 0.5,
  #                 nudge_x = .5,
  #                 nudge_y = c(.1, -.05)
  # ) +
  facet_wrap(.~str_to_title(domain)) +
  theme_minimal() +
  labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
       title = "Behaviour") +
  theme(strip.placement = "outside",
        legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
        legend.margin = margin(-.5,0,0,0, unit="cm"),
        legend.spacing.y = unit(0.15, 'cm'),
        legend.key.width = unit(1, "cm"),
        legend.key.size = unit(.3, "cm"),
        # plot.tag.position = c(0,.8),
        legend.text = element_text(size = 10, color = "grey20"),
        text = element_text(size = 12, color = "grey40"),
        axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
        axis.text.x = element_text( hjust=c(0,1)),
        title = element_text(size = 12, color = "grey20"),
        plot.tag  = element_text( size = 12, face = "bold", color = "grey20"),
        panel.spacing = unit(.5, "lines"),
        plot.title = element_blank(),
        panel.grid = element_blank(),
        plot.title.position = "plot",
        plot.margin = margin(b = 5, r = 5, l = 5),
        panel.background = element_rect(color = "grey75", fill = NA, size = .4),
        strip.text = element_text(face = "bold"),
        strip.text.y = element_text(size = 12),
        strip.text.x = element_text(size = 12)) +
  scale_fill_manual(values = c("#502a7b","#502a7c" , "#502a7a" )) +
  guides(color = guide_legend(override.aes = list(size = .75)),
         fill =  guide_legend(override.aes = list(size = .75)),
         size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
  scale_y_continuous(breaks = seq(0,1,0.25)) +
  scale_x_continuous(breaks = c(0,1,2,3,4,5))

2.2.4 Plot correlation and interval by item

nd <- crossing(domain= NA,  
               sei = 0.1,
               item= unique(df$item), 
               event= NA,
               formulation= NA,
               age_dec_c= 0,
               age_dec_c2= 0,
               female_c= 0, 
               interval_val=seq(0,5, by = .1))


epred_draws_df <- nd %>% 
  add_epred_draws(fit_masc_m3, re_formula = NA)

# epred_draws_dom <- epred_draws_df %>%
#   group_by(interval_val, health_subdomain) %>%
#   mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
#   pivot_wider(names_from = .width, values_from = c(.lower,.upper))

epred_draws_df  <- epred_draws_df  %>%
  mutate(item = str_to_title(as.character(item)))


ggplot(epred_draws_df) +
  stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) + 
  geom_point(data= df, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
  # geom_line(data = epred_draws_agg,
  #           aes(x = time_diff_dec*10, y = .epred),
  #           color = "grey95",
  #           size = .5) +
  # geom_line(data = epred_draws_dom, 
  #           aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
  #           color = "#e07f00",
  #           linewidth = .25) +
  # geom_text_repel(data = lbl_dot_df, 
  #                 aes(x = time_diff_dec*10, y = .epred, label = label),
  #                 family = "Source Sans 3", size = 2.5,
  #                 min.segment.length = 0,
  #                 segment.color = "grey50",
  #                 segment.size = .25,
  #                 box.padding = 0.5,
  #                 nudge_x = .5,
  #                 nudge_y = c(.1, -.05)
  # ) +
  facet_wrap(.~str_to_title(item)) +
  theme_minimal() +
  labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
       title = "Behaviour") +
  theme(strip.placement = "outside",
        legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
        legend.margin = margin(-.5,0,0,0, unit="cm"),
        legend.spacing.y = unit(0.15, 'cm'),
        legend.key.width = unit(1, "cm"),
        legend.key.size = unit(.3, "cm"),
        # plot.tag.position = c(0,.8),
        legend.text = element_text(size = 10, color = "grey20"),
        text = element_text(size = 12, color = "grey40"),
        axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
        axis.text.x = element_text( hjust=c(0,1)),
        title = element_text(size = 12, color = "grey20"),
        plot.tag  = element_text( size = 12, face = "bold", color = "grey20"),
        panel.spacing = unit(.5, "lines"),
        plot.title = element_blank(),
        panel.grid = element_blank(),
        plot.title.position = "plot",
        plot.margin = margin(b = 5, r = 5, l = 5),
        panel.background = element_rect(color = "grey75", fill = NA, size = .4) ,
        strip.text = element_text(face = "bold"),
        strip.text.y = element_text(size = 12),
        strip.text.x = element_text(size = 12))+ 
  scale_fill_manual(values = c("#502a7b","#502a7c" , "#502a7a" )) +
  guides(color = guide_legend(override.aes = list(size = .75)),
         fill =  guide_legend(override.aes = list(size = .75)),
         size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0,5)) +
  scale_y_continuous(breaks = seq(0,1,0.25)) +
  scale_x_continuous(breaks = c(0,1,2,3,4,5))

2.2.5 Plot correlation and interval by formulation

nd <- crossing(domain= NA,  
               sei = 0.1,
               item= NA, 
               event= NA,
               formulation= unique(df$formulation),
               age_dec_c= 0,
               age_dec_c2= 0,
               female_c= 0, 
               interval_val=seq(0,5, by = .1))


epred_draws_df <- nd %>% 
  add_epred_draws(fit_masc_m3, re_formula = NA)

# epred_draws_dom <- epred_draws_df %>%
#   group_by(interval_val, health_subdomain) %>%
#   mean_hdci(.epred,.width = c(.95,.8,.5)) %>%
#   pivot_wider(names_from = .width, values_from = c(.lower,.upper))


# Set order and capitalise first letter
epred_draws_df$formulation <- str_to_title(epred_draws_df$formulation)

epred_draws_df$formulation<- factor(epred_draws_df$formulation, levels = c("Risk","Concern", "Worry" ,"Threat"))

df_plot <- df %>% 
  mutate(formulation = str_to_title(formulation)) 

df_plot$formulation<- factor(df_plot$formulation, levels = c("Risk","Concern", "Worry" ,"Threat"))

# Plot
ggplot(epred_draws_df) +
  stat_lineribbon(alpha = 1/4, point_interval = "mean_hdci", aes(x = interval_val, y = .epred)) + 
  geom_point(data= df_plot, aes(x=interval_val, y= cor_val), size = 1.5, color = "grey20", alpha = 0.7) +
  # geom_line(data = epred_draws_agg,
  #           aes(x = time_diff_dec*10, y = .epred),
  #           color = "grey95",
  #           size = .5) +
  # geom_line(data = epred_draws_dom, 
  #           aes(x = time_diff_dec*10, y = .epred, linetype = domain_name),
  #           color = "#e07f00",
  #           linewidth = .25) +
  # geom_text_repel(data = lbl_dot_df, 
  #                 aes(x = time_diff_dec*10, y = .epred, label = label),
  #                 family = "Source Sans 3", size = 2.5,
  #                 min.segment.length = 0,
  #                 segment.color = "grey50",
  #                 segment.size = .25,
  #                 box.padding = 0.5,
  #                 nudge_x = .5,
  #                 nudge_y = c(.1, -.05)
  # ) +
  facet_wrap(.~formulation) +
  theme_minimal() +
  labs(y = "Retest Correlation", x = "Retest Interval (Years)", color = "", linetype = "", fill = "",
       title = "Behaviour") +
  theme(strip.placement = "outside",
        legend.position = "none", # c(0,0) bottom left, c(1,1) top-right.
        legend.margin = margin(-.5,0,0,0, unit="cm"),
        legend.spacing.y = unit(0.15, 'cm'),
        legend.key.width = unit(1, "cm"),
        legend.key.size = unit(.3, "cm"),
        # plot.tag.position = c(0,.8),
        legend.text = element_text(size = 10, color = "grey20"),
        text = element_text(size = 12, color = "grey40"),
        axis.text.y = element_text( vjust=seq(0,1, length.out = 5)),
        axis.text.x = element_text( hjust=c(0,1)),
        title = element_text(size = 12, color = "grey20"),
        plot.tag  = element_text( size = 12, face = "bold", color = "grey20"),
        panel.spacing = unit(.5, "lines"),
        plot.title = element_blank(),
        panel.grid = element_blank(),
        plot.title.position = "plot",
        plot.margin = margin(b = 5, r = 5, l = 5),
        panel.background = element_rect(color = "grey75", fill = NA, size = .4),
        strip.text = element_text(face = "bold"),
        strip.text.y = element_text(size = 12),
        strip.text.x = element_text(size = 12)) +
  scale_fill_manual(values = c("#502a7b","#502a7c","#502a7a" ))+
  guides(color = guide_legend(override.aes = list(size = .75)),
         fill =  guide_legend(override.aes = list(size = .75)),
         size = "none", linetype = guide_legend(override.aes = list(size = .75))) +
  coord_cartesian(ylim = c(0, 1), xlim = c(0,5))+
  scale_y_continuous(breaks = seq(0,1,0.25)) +
  scale_x_continuous(breaks = c(0,1,2,3,4,5))

2.3 Prediction Parameters

2.3.1 By Domain

pred_df_domain <- NULL

for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= unique(df$domain),  
                 sei = 0.1,
                 item=NA, 
                 event=NA,
                 formulation=NA,
                 age_dec_c=0,
                 age_dec_c2=0,
                 female_c= 0, 
                 interval_val=0)
  
  
  
  
  fit_nlpar_domain <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA)    
  

 
  fit_nlpar_domain <- fit_nlpar_domain %>%
    group_by(domain) %>% 
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "domain",
           x = domain)%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  pred_df <- fit_nlpar_domain
  

  pred_df_domain <- bind_rows(pred_df, pred_df_domain) 
}



pred_df_domain <- pred_df_domain %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.2 By Item

pred_df_item <- NULL

  
  
  for (curr_nlpar in c("stabch","rel","change")) {
    
    
    nd <- crossing(domain= NA,  
                   sei = 0.1,
                   item= unique(df$item), 
                   event=NA,
                   formulation=NA,
                   age_dec_c=0,
                   age_dec_c2=0,
                   female_c= 0, 
                   interval_val=0)
    
    
    fit_nlpar_item <- nd %>% 
      add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
    

  
    
    
    fit_nlpar_item <- fit_nlpar_item %>%
      group_by(item) %>% 
      mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
      pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
      mutate(nlpar = curr_nlpar,
             estimate = .epred,
             categ = "item",
             x = item)%>% 
      select(categ, x, nlpar, estimate, dplyr::contains("er_"))
    
    
    
    pred_df <- fit_nlpar_item
    
    
    pred_df_item <- bind_rows(pred_df, pred_df_item) 
  }
  

pred_df_item <- pred_df_item %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.3 By Event

pred_df_event <- NULL



for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= NA,  
                 sei = 0.1,
                 item= NA, 
                 event=unique(df$event),
                 formulation=NA,
                 age_dec_c=0,
                 age_dec_c2=0,
                 female_c= 0, 
                 interval_val=0)
  
  
  fit_nlpar_event <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
  
  
  
  
  
  fit_nlpar_event <- fit_nlpar_event %>%
    group_by(event) %>% 
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "event",
           x = event)%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  
  pred_df <- fit_nlpar_event
  
  
  pred_df_event <- bind_rows(pred_df, pred_df_event) 
}

pred_df_event <- pred_df_event %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.4 By Formulation

pred_df_formulation <- NULL



for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= NA,  
                 sei = 0.1,
                 item= NA, 
                 event=NA,
                 formulation= unique(df$formulation),
                 age_dec_c=0,
                 age_dec_c2=0,
                 female_c= 0, 
                 interval_val=0)
  
  
  fit_nlpar_formulation <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
  
  
  
  
  
  fit_nlpar_formulation <- fit_nlpar_formulation %>%
    group_by(formulation) %>% 
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "formulation",
           x = formulation)%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  
  pred_df <- fit_nlpar_formulation
  
  
  pred_df_formulation <- bind_rows(pred_df, pred_df_formulation) 
}

pred_df_formulation <- pred_df_formulation %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.5 By Age

pred_df_age <- NULL



for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= NA,  
                 sei = 0.1,
                 item= NA, 
                 event=NA,
                 formulation=NA,
                 age_dec_c= c(-2,0,2,4),
                 age_dec_c2= c(4, 0, 4, 16),
                 female_c= 0, 
                 interval_val=0)
  
  
  fit_nlpar_age <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
  
  
  
  
  
  fit_nlpar_age <- fit_nlpar_age %>%
    group_by(age_dec_c) %>% 
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "age",
           x = case_when(age_dec_c== 4 ~ "80+", age_dec_c== 2 ~ "60-79", age_dec_c== 0 ~ "40-59", age_dec_c== -2 ~ "20-39"))%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  
  pred_df <- fit_nlpar_age
  
  
  pred_df_age <- bind_rows(pred_df, pred_df_age) 
}

pred_df_age <- pred_df_age %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.6 By Gender

pred_df_female <- NULL



for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= NA,  
                 sei = 0.1,
                 item= NA, 
                 event=NA,
                 formulation=NA,
                 age_dec_c= 0,
                 age_dec_c2= 0,
                 female_c= c(-0.5, 0.5), 
                 interval_val=0)
  
  
  fit_nlpar_female <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
  
  
  
  
  
  fit_nlpar_female <- fit_nlpar_female %>%
    group_by(female_c) %>% 
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "gender",
           x = ifelse(female_c== -0.5, "male", "female"))%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  
  pred_df <- fit_nlpar_female
  
  
  pred_df_female <- bind_rows(pred_df, pred_df_female) 
}

pred_df_female <- pred_df_female %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))

2.3.7 Overall

pred_df_overall <- NULL



for (curr_nlpar in c("stabch","rel","change")) {
  
  
  nd <- crossing(domain= NA,  
                 sei = 0.1,
                 item= NA, 
                 event=NA,
                 formulation=NA,
                 age_dec_c= 0,
                 age_dec_c2= 0,
                 female_c= 0, 
                 interval_val=0)
  
  
  fit_nlpar_overall <- nd %>% 
    add_epred_draws(fit_masc_m3, nlpar = curr_nlpar, re_formula = NA) 
  
  
  
  
  
  fit_nlpar_overall <- fit_nlpar_overall %>%
    mean_hdci(.epred,.width = c(.95,.8,.5)) %>% 
    pivot_wider(names_from = .width, values_from = c(.lower,.upper)) %>% 
    mutate(nlpar = curr_nlpar,
           estimate = .epred,
           categ = "all",
           x = "overall")%>% 
    select(categ, x, nlpar, estimate, dplyr::contains("er_"))
  
  
  
  pred_df <- fit_nlpar_overall
  
  
  pred_df_overall <- bind_rows(pred_df, pred_df_overall) 
}


pred_df_overall <- pred_df_overall %>% 
  mutate(nlpar = case_when(nlpar == "rel" ~ "Reliability",
                           nlpar == "change" ~ "Change",
                           nlpar == "stabch" ~"Stab. Change"))


pred_df <- bind_rows(pred_df_domain, pred_df_item, pred_df_event, pred_df_formulation ,pred_df_age, pred_df_female, pred_df_overall)

2.4 Parameters Table

library(gt)

# Selecting only the specified columns from pred_df
pred_df_selected <- pred_df[, c("categ", "x", "nlpar", "estimate", ".lower_0.95", ".upper_0.95")]

# Create a gt table
table_gt <- gt(data = pred_df_selected,
               rowname_col = "x",
               caption = "Prediction Results",
               rownames_to_stub = TRUE)

# Formatting
table_gt <- table_gt %>%
  tab_header(title = "Estimates and Confidence Intervals") %>%
  fmt_number(columns = c("estimate", ".lower_0.95", ".upper_0.95"), decimals = 3) %>%
  cols_align(align = "center") 

table_gt
Prediction Results
Estimates and Confidence Intervals
categ x nlpar estimate .lower_0.95 .upper_0.95
1 domain health Change 0.459 0.069 0.881
2 domain other Change 0.428 0.039 0.850
3 domain health Reliability 0.775 0.657 0.897
4 domain other Reliability 0.772 0.648 0.900
5 domain health Stab. Change 0.480 0.090 0.898
6 domain other Stab. Change 0.539 0.111 0.945
7 item multiple Change 0.425 0.054 0.872
8 item single Change 0.459 0.038 0.914
9 item multiple Reliability 0.745 0.624 0.864
10 item single Reliability 0.795 0.640 0.942
11 item multiple Stab. Change 0.408 0.033 0.862
12 item single Stab. Change 0.609 0.161 0.973
13 event none Change 0.394 0.042 0.826
14 event event Change 0.490 0.094 0.923
15 event none Reliability 0.838 0.739 0.932
16 event event Reliability 0.694 0.553 0.848
17 event none Stab. Change 0.461 0.055 0.886
18 event event Stab. Change 0.559 0.153 0.950
19 formulation risk Change 0.269 0.001 0.787
20 formulation concern Change 0.492 0.078 0.930
21 formulation threat Change 0.555 0.127 0.943
22 formulation worry Change 0.504 0.103 0.907
23 formulation risk Reliability 0.656 0.541 0.784
24 formulation concern Reliability 0.773 0.542 0.989
25 formulation threat Reliability 0.711 0.545 0.875
26 formulation worry Reliability 0.879 0.784 0.970
27 formulation risk Stab. Change 0.730 0.203 1.000
28 formulation concern Stab. Change 0.467 0.053 0.893
29 formulation threat Stab. Change 0.396 0.024 0.816
30 formulation worry Stab. Change 0.407 0.038 0.830
31 age 20-39 Change 0.583 0.011 1.000
32 age 40-59 Change 0.651 0.034 1.000
33 age 60-79 Change 0.710 0.029 1.000
34 age 80+ Change 0.737 0.011 1.000
35 age 20-39 Reliability 0.744 0.523 0.983
36 age 40-59 Reliability 0.766 0.574 1.000
37 age 60-79 Reliability 0.781 0.576 1.000
38 age 80+ Reliability 0.790 0.556 1.000
39 age 20-39 Stab. Change 0.496 0.000 0.998
40 age 40-59 Stab. Change 0.558 0.000 0.998
41 age 60-79 Stab. Change 0.614 0.000 0.999
42 age 80+ Stab. Change 0.646 0.000 1.000
43 gender male Change 0.435 0.083 0.830
44 gender female Change 0.448 0.093 0.835
45 gender male Reliability 0.752 0.603 0.893
46 gender female Reliability 0.793 0.672 0.911
47 gender male Stab. Change 0.523 0.141 0.892
48 gender female Stab. Change 0.496 0.103 0.857
49 all overall Change 0.439 0.114 0.805
50 all overall Reliability 0.775 0.667 0.894
51 all overall Stab. Change 0.510 0.159 0.861

2.5 Parameters Plot

pred_df <- pred_df %>%
  mutate(categ = str_to_title(as.character(categ)),
         x = str_to_title(as.character(x)))


pred_df$nlpar <- factor(pred_df$nlpar, levels = c("Reliability", "Change", "Stab. Change"))


pred_df$x <- factor(pred_df$x, levels = c("Overall", "None","Event", "Other", "Health",  "Single", "Multiple", "Threat", "Concern", "Worry", "Risk", "80+","60-79", "40-59", "20-39", "Male", "Female"))

pred_df$categ <- factor(pred_df$categ, levels = c("All", "Event", "Domain",  "Item", "Formulation"  ,"Age", "Gender"))


p_nlpar <- pred_df %>% ggplot() +
  geom_crossbar(aes(xmin = .lower_0.95, x = estimate, 
                    xmax = .upper_0.95, y = x),
                fill = "white", color = "NA",
                linewidth = .15,width = 0.25, alpha =  1) +
  geom_crossbar(aes(xmin = .lower_0.95, x = estimate, 
                    xmax = .upper_0.95, y = x),
                fill = "#502a7a", color = "NA",
                linewidth = .15,width = 0.25, alpha =  .3) +
  geom_crossbar(aes(xmin = .lower_0.8, x = estimate,
                    xmax = .upper_0.8, y = x),
                fill = "white", color = "NA",
                linewidth = .15,width = 0.25, alpha =  1) +
  geom_crossbar(aes(xmin = .lower_0.8, x = estimate,
                    xmax = .upper_0.8, y = x),
                fill = "#502a7a",color = "NA",
                linewidth = .15,width = 0.25, alpha =  .6) +
  geom_crossbar(aes(xmin = .lower_0.5, x = estimate, 
                    xmax = .upper_0.5, y = x),
                fill = "white",color = "NA",
                linewidth = .15,width = 0.25, alpha =  1) +
  geom_crossbar(aes(xmin = .lower_0.5, x = estimate, 
                    xmax = .upper_0.5, y = x),
                fill = "#502a7a",color = "NA",
                linewidth = .15,width = 0.25, alpha =  .9) +
  geom_point(aes(x = estimate, y =x),
             fill = "white", color = "grey20",
             shape = 21, 
             stroke = .25, 
             size = 2) +
  facet_grid(categ ~ nlpar, scales = "free_y", space = "free", switch = "y") + 
  # scale_x_continuous(limits = c(0,1), expand = c(0,0)) +
  scale_x_continuous(limits = c(0, 1), expand = c(0, 0), breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1")) +
  scale_y_discrete(position = "left") +
  theme_minimal() +
  geom_rect(data = subset(pred_df, x %in% c("Overall")), 
            fill = NA, colour = "black", size = .75, xmin = -Inf,xmax = Inf,
            ymin = -Inf,ymax = Inf) +
  theme(panel.grid = element_blank(),
        legend.position = "none",
        plot.title.position = "plot",
        strip.placement = "outside",
        #strip.text.y = element_blank(),
        axis.title.x = element_text(size = 10, color = "grey20"),
        plot.margin = margin(b = 10, t = 10, r = 15, l = 0),
        panel.spacing.x = unit(.3, "cm"),
        # plot.tag.position = c(0.025,.9),
        # panel.grid.major.x = element_line(linewidth = .2, color = "grey80", linetype = "solid"),
        panel.background = element_rect(linewidth = .25, color = "grey50", fill = "NA"),
        #plot.background = element_rect(linewidth = .25, color = "grey40", fill = "NA"),
        strip.text =  element_text( face = "bold"),
        strip.text.y = element_text(size = 10),
        strip.text.x = element_text(size = 10),
        plot.tag  = element_text( size = 10, face = "bold", color = "grey20"),
        plot.title = element_blank(),
        title = element_text( face = "bold", size = 10),
        axis.text.x =  element_markdown( color = "black", size = 10, hjust=c(0,.5, 1)),
        axis.text.y.left =  element_markdown( angle = 0, hjust = 1, color = "grey20", size = 10), # hjust = c(0,.5,.5,.5,1)
        text = element_text()) +
  labs(y = "", x = "Parameter Estimate", title = "Propensity") 


p_nlpar

ggsave(plot = p_nlpar,
       filename = here::here("analysis", "masc_pred.jpeg"), 
       dpi = 300, width = 28, height = 30, units = "cm")